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Abstract. I give an introduction to the basic concepts of fluid dynamics as applied 
to the dynamical description of relativistic nuclear collisions. 



1 Introduction and Conclusions 

Modelling the dynamic evolution of nuclear collisions in terms of fluid dynamics 
has a long-standing tradition in heavy- ion physics, for a review see [1, 2, 3]. One of 
the main reasons is that one essentially does not need more information to solve 
the equations of motion of ideal fluid dynamics than the equilibrium equation of 
state of matter under consideration. Once the equation of state is known (and 
an initial condition is specified), the equations of motion uniquely determine 
the dynamics of the collision. Knowledge about microscopic reaction processes 
is not required. This becomes especially important when one wants to study 
the transition from hadron to quark and gluon degrees of freedom, as predicted 
by lattice simulations of quantum chromodynamics (QCD) [4]. The complicated 
deconfinement or hadronization processes need not be known in microscopic 
detail, all that is necessary is the thermodynamic equation of state as computed 
by e.g. lattice QCD. This fact has renewed interest in fluid dynamics to study the 
effects of the deconfinement and chiral symmetry restoration transition on the 
dynamics of relativistic nuclear collisions. Such collisions are presently under 
intense experimental investigation at CERN's SPS and Brookhaven National 
Laboratory's AGS and (beginning Fall 1999) RHIC accelerators. 

In this set of lectures I give an overview over the basic concepts and notions 
of relativistic fluid dynamics as applied to the physics of heavy-ion collisions. 
The aim is not to present a detailed and complete review of this field, but to 
provide a foundation to understand the literature on current research activities 
in this field. This has the consequence that the list of references is far from com- 
plete, that I will not make any attempt to compare to actual experimental data, 
and that some interesting, but more applied topics (such as transverse collective 
flow) will not be discussed here. In Section 2, I discuss the basic concepts of 
relativistic fluid dynamics. First, I present a derivation of the fluid-dynamical 
equations of motion. A priori, there are more unknown functions than there are 
equations, and one has to devise approximation schemes in order to close the 
set of equations of motion. The most simple is the ideal fluid approximation. 
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which simply discards some of the unknown functions. Another one is the as- 
sumption of small deviations from local thermodynamical equilibrium, which 
leads to the equations of dissipative fluid dynamics. A brief discussion of multi- 
fluid models concludes this section. In Section 3 I discuss numerical aspects of 
solution schemes for ideal relativistic fluid dynamics. Section 4 is devoted to a 
discussion of one-dimensional solutions of ideal fluid dynamics. After present- 
ing a classification of possible wave patterns in one spatial dimension, for both 
thermodynamicall normal as well as anomalous matter, I discuss the expansion 
of semi-inflnite matter into vacuum. This naturally leads to the discussion of 
the Landau model for the one- dimensional expansion of a finite slab of matter. 
The Landau model is historically the first fluid-dynamical model for relativistic 
nuclear collisions. However, more realistic is, at least for ultrarelativistic collision 
energies, the so-called the Bjorken model which is subsequently prc;sented. The 
main result of this section is the time delay in the expansion of the system due 
to the softening of the equation of state in a phase transition region. This may 
have potential experimental consequences for nuclear collisions at RHIC energies, 
where one wants to study the transition from hadron to quark and gluon degrees 
of freedom. Finally, Section 5 concludes this set of lectures with a discussion on 
how to decouple particles from the fluid evolution in the so-called "freeze-out" 
process and compute experimentally observable quantities like single inclusive 
particle spectra. 

Units are h ~ c = = I. The metric tensor is g^'^ = diag(+, — , — . — ). 
Upper greek indices are contravariant, lower greek indices covariant. The scalar 
product of two 4- vectors , is denoted by g^i, V = a'^h^ = a ■ h. The 
latter notation is also used for the scalar product of two 3-vectors a, b, a • b. 

2 The Basics 

In this section, I flrst derive the conservation equations of relativistic fluid dy- 
namics. If there are n conserved charges in the fluid, there are 4 -|- n conser- 
vation equations: 1 for the conservation of energy, 3 for the conservation of 
3-momentum, and n for the conservation of the respective charges. In the gen- 
eral case, however, there are 10-|-4n independent variables: the 10 independent 
components of the energy-momentum tensor (which is a symmetric tensor of 
rank 2), and the 4 independent components of the 4- vectors of the n charge 
currents. Thus, the system of fluid-dynamical equations is not closed, and one 
requires an approximation in order to solve it. 

The simplest approximation is the ideal fluid assumption which reduces the 
number of unknown functions to 5 -|- n. The equation of state of the fluid then 
provides the final equation to close the system of conservation equations and to 
solve it uniquely. Another approximation is the assumption of small deviations 
from an ideal fluid and leads to the equations of dissipative fluid dynamics. In 
this approximation one provides an additional set of 6-1-3 n equations to close the 
set of equations of motion. Finally, I also briefly discuss multi-fluid-dynamical 
models. 



Fluid Dynamics for Relativistic Nuclear Collisions 



3 



2.1 The Conservation Equations 

Fluid dynamics is equivalent to the conservation of energy, momentum, and net 
charge number. Consider a single fluid characterized locally in space-time by its 
energy- momentum tensor T^^"(x) and by the n conserved net charge currents 
N-^{x) , i = l,...,n. (Conserved charges are for example the electric charge, 
baryon number, strangeness, charm, etc.) Consider now an arbitrary hyper- 
surface E in 4-dimensional space-time. The tangent 4-vector on this surface is 
S'' (.'/;) . The normal vector on a surface element dS of E is denoted by dZ'^ (x) . 
By definition, dU ■ S = 0. The amount of net charge of type i and of energy and 
momentum flowing through the surface element dS is given by 

dNi = dIJ ■ Ni , i = l,...,n , (1) 
dP^" = dE^T'"' , i^ = 0,...,3 . (2) 

Now consider an arbitrary space-time volume V4 with a closed surface S. If there 
are neither sources nor sinks of net charge and energy-momentum inside S, one 
has 

(j) dE-Ni = , i = l,...,n , (3) 
^ dZ"^ T**^ = , At = 0, . . . , 3 . (4) 

Gauss theorem then leads immediately to the global conservation of net charge 
and energy-momentum: 

[ d^xa^TVf^O, i = l,...,n , (5) 
[ d'^xd^T^"' = , ij. = 0,...,3 . (6) 

For arbitrary V4, however, one has to require that the integrands in (5,6) vanish, 
which leads to local conservation of net charge and energy-momentum: 

d^Nt = 0, t=l,...,n , (7) 
d^T^'' = 0, iy = 0,...,3 . (8) 

These are the equations of motion of relativistic fluid dynamics [5]. Note that 
there are A-\-n equations, but 10+4 n independent unknown functions T^^(x) , N^{x). 
{T^^ is a symmetric tensor of rank 2 and therefore has 10 independent compo- 
nents, the are 4-vectors with 4 independent components.) Therefore, the 
system of fluid-dynamical equations is a priori not closed and cannot be solved 
in complete generality. One requires additional assumptions to close the set of 
equations. One possibility is to reduce the number of unknown functions, an- 
other is to provide 6 -|- 3 n additional equations of motion which determine all 
unknown functions uniquely. Both possibilities will be discussed in the following 
subsections. 
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2.2 Tensor Decomposition and Choice of Frame 

Before discussing approximations to close the system of conservation equations, 
it is convenient to perform a tensor decomposition of , T''^ with respect to 
an arbitrary, time-like, normalized 4- vector u'^ , u - u = 1. The projector onto the 
3-space orthogonal to u'* is denoted by 

^^,u ^ _ ^My^ ^ = , Af^^A"^ = A^"" . (9) 

Then the tensor decomposition reads: 

A^f = n, + i^^ , (10) 

T^"" = eu''u'' -pA'''' +q''u'' + q''u'' + TT'''' , (11) 

where 

n, = N,-u (12) 

is the net density of charge of type i in the frame where = (1,0) (subsequently 
denoted as the local rest frame, LRF), 

= (13) 

is the net flow of charge of type i in the LRF, 

e = u^T'^'u, (14) 

is the energy density in the LRF, 

(15) 

is the isotropic pressure in the LRF, 

< = A^'°'Tc,0u'^ (16) 
is the flow of energy or heat flow in the LRF, and 

.a/3 (^7) 



^^a^^a':,)-\a^'^a^^ 



is the stress tensor in the LRF. Note that the particular projection (17) is trace- 
free. (The trace of the projection A'^T'^^A'^ is absorbed in the definition of p.) 
The tensor decomposition replaces the original 10 + 4n unknown functions by 
an equal number of new unknown functions rij (n), (3n), e(l), p(l), (3), 
and ttI^" (5). 

So far, u^^ is arbitrary. However, one can give it a physical meaning by choos- 
ing it either to be 
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or (which is an implicit definition) 



< ^ , '^^^ „ . (19) 



is the physical 4- velocity of the flow of net charge 
i. The LRF is then the local rest frame of the flow of net charge i, i.e., the frame 
where A^^^ = {N^, 0). In this frame, there is obviously no flow of charge i, = 0, 
and = m. This LRF is called Eckart frame. Note, however, that not all net 
charges need to flow with the same velocity, might be ^ for j ^ i. The 
number of unknown functions is still 10 + 4n, since the 3 previously unknown 
functions t'^ have been merely replaced by the 3 independent components of 
(ue ■ Ue = 1!), which now have to be determined dynamically from the 
conservation equation for Nj^ . 

The second choice means that is the physical 4- velocity of the energy flow. 
The LRF is the local rest frame of the energy flow. It is obvious that in this frame 
q'^ = 0. This frame is called Landau fram,e. The number of unknown functions 
is still 10 + 4n, since the 3 previously unknown functions q'^ have been merely 
replaced by the 3 independent components of ('"l • = 1!); which now have 
to be determined dynamically from the conservation equation for T^^'^. Other 
choices of rest frames are also possible, for a discussion, see [6]. 



2.3 Ideal Fluid Dynamics 

Consider an ideal gas in local thermodynamical equilibrium. The single-particle 
phase space distribution for fermions or bosons then reads 

/o(^' ^) = ei^pik-u{x)- fi{x))/T{x)±l ' 

where u^{x) is the local average 4-velocity of the particles, /i(x) and T{x) are 
local chemical potential and temperature, and g counts internal degrees of free- 
dom (spin, isospin, color, etc.) of the particles. The chemical potential of the 
particles is defined as /i = X]r=i where /ij are the chemical potentials which 
control the net number of charge of type i, and qi is the individual charge of 
type i carried by a particle. The chemical potential for antiparticles is /x = — /z 
(in thermodynamical equilibrium) . Let us define the single-particle phase space 
distribution for antiparticles by fa{fl) = /o(— At)- 

The kinetic definitions of the net current of charge of type i and of the 
energy-momentum tensor are [6] 

N^{x) ^q^J^k'^ [fo{k, x) - fo{k, x)] , (21) 

/d^k 
— fc^r [fo{k,x) + fo{k,x)] , (22) 
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where E = Vk^ + is the on-shell energy of the particles and m their rest 
mass. Inserting (20) one computes 




where 




is the thermodynamic net number density of charge of type i of an ideal gas, 
and the Fermi-Dirac or Bose-Einstein distribution was denoted by n{E) = 



l/{exp[{E - n)/T] ± 1), n{E) = l/(exp[(£; + /x)/T] ± 1). Furthermore, 



is the thermodynamic ideal gas pressure. The form (23,24) impHes that for an 
ideal gas in local thermodynamical equilibrium the functions vf = = tt^^ = 0, 
i.e., there is no flow of charge or heat with respect to the particle flow velocity 
u'*, and there are no stress forces. This implies furthermore (and can be con- 
firmed by an explicit calculation) that for an ideal gas in local thermodynamical 
equilibrium u^ = u'^ = u^, i.e., Eckart's and Landau's choice of frame coincide 
with the local rest frame of particle flow. 

This consideration of an ideal gas in local thermodynamical equilibrium 
serves as a motivation for the so-called ideal fluid approximation. In this ap- 
proximation, one starts on the macroscopic level of fluid variables iVf , Ti^" and 
a priori takes them to be of the form (23) and (24). The corresponding fluid is 
referred to as an ideal fluid. Without any further assumption, however, the corre- 
sponding system of 4 + n equations of motion contains 5 -I- n unknown functions, 
e,p, u^, and Ui, i = 1, . . . ,n. One therefore has to specify an equation of state for 
the fluid, for instance (and most commonly taken) of the form p(e, ni, . . . ,n„). 
This closes the system of equations of motion. 

The equation of state is the only place where information enters about the na- 
ture of the particles in the fluid and the microscopic interactions between them. 
Usually, the equation of state for the fluid is taken to be the thermodynamic 
equation of state, as computed for a system in thermodynamical equilibrium. 
The process of closing the system of equations of motions by assuming a ther- 
modynamic equation of state therefore involves the implicit assumption that the 
fluid is in local thermodynamical equilibrium. It is important to note, however, 
that the explicit form of the equation of state is completely unrestricted, for 
instance it can have anomalies like phase transitions. 





is the thermodynamic ideal gas energy density, and 




(27) 
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The ideal fluid approximation therefore allows to consider a wider class of 
systems than just an ideal gas in local thermodynamical equilibrium, which 
served as a motivation for this approximation. An ideal gas has a very specific 
equation of state without any anomalies and is given by (27) which defines 
p{T, Hi, . . . , fXn) (which in turn allows to determine all other thermodynamic 
functions from the first law and the fundamental relation of thermodynamics, 
and thus to specify p(e, ni, . . . , n„), see the following remarks). 

I close this subsection with three remarks. The first concerns the notion of 
an equation of state which is com,plete in the thermodynamic sense. Such an 
equation of state allows (by definition) to determine, for given values of the 
independent thermodynamic variables, all other thermodynamic functions from 
the first law of thermodynamics (or one of its Legendre transforms) 

1 " 

d.= ide-^^n, , (28) 

i=l 

s being the entropy density, and from the fundamental relation of thermody- 
namics 

n 

e+p = Ts + Y,l^ini . (29) 

Obviously, for independent thermodynamic variables e, ni, . . . , n„, an equation 
of state of the form s(e, ni, . . . , n„) is complete in this sense, since partial diff'er- 

entiation of this fimction yields, from (28), the fimctions 1/T, /ii/T, . . ., i-in/T. 
Then, the fundamental relation (29) yields the last unknown thermodynamic 
function, p. 

Another example of a complete equation of state is p(T, /xi, . . . , /i„), since 
the (multiple) Legendre transform of (28) reads 

n 

dp = s dT + ^ nj d/ii (30) 

i=l 

(which is also known as the Gibbs-Duhem relation), such that the thermody- 
namic functions s, m, . . . , n„ can be determined from partial differentiation of 
p. The last unknown thermodynamic function, e, can then be determined from 
(29). 

The equation of state p(e, ni, . . . , n„) is, however, not a complete equation of 
state in the thermodynamic sense. Partial differentiation of this function yields 
thermodynamic functions dp/de, dp/drii, i = 1, . . . ,n, which in general do not 
allow to infer the values of T, s, and fii, i = 1, . . . ,n. 

The second remark concerns the assumption of local thermodynamical equi- 
librium. In order to achieve local thermodynamical equilibrium, spatio-temporal 
variations of the macroscopic fluid flelds have to be small compared to micro- 
scopic reaction rates which drive the system (locally) towards thermodynami- 
cal equilibrium. A quantity that characterizes spatio-temporal variations of the 
macroscopic fields is the so-called expansion scalar 6 = d ■ u. It determines the 
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(local) rate of expansion of the fluid. Microscopic reaction rates arc essentially 
given by the product of cross section and local particle density, F ~ crn. The 
criterion for local thermodynamical equilibrium then reads 

r>6', or cr>6'/n . (31) 

The third remark concerns entropy production. In ideal fluid dynamics, the 
entropy current is defined as 

S^ = su^ . (32) 

Taking the projection of energy-momentum conservation in the direction of Ui, 
one derives 

= d^T^"" = e+{e + p)e , (33) 

where d = u ■ da is a. comoving time derivative and where use has been made 

of the fact that is normalized, i.e., 9^(m • u) = 0. With the first law of 
thermodynamics (28) and the fundamental relation of thermodynamics (29) one 
rewrites this as 

n 

T{s + se)+J2iJ.i{ni+nie)=0 . (34) 

i=l 

Finally, employing net charge conservation d ■ Ni = hi + riiO = yields 

s + se = d-S = , (35) 

i.e., the entropy current is conserved in ideal fiuid dynamics. As we shall sec in 
one of the following section, however, this proof only holds where the partial 
derivatives in these equations are well-defined, i.e., for continuous solutions of 
ideal fluid dynamics. Discontinuous solutions will in fact be shown to produce 
entropy. 



2.4 Dissipative Fluid Dynamics 

In dissipative fluid dynamics one docs not set i^f , g^, tt^" a priori to zero, but 
specifies them through additional equations. There are two ways to obtain the 
latter. The first is phenomenological and starts from the second law of thermo- 
dynamics, i.e., the principle of non-decreasing entropy, 

d-S>0 . (36) 

The second way resorts to kinetic theory to derive the respective equations. In 
principle, both ways require the additional assumption that deviations from local 
thermodynamical equilibrium are smaU. To make this statement more concise, 
let us introduce the equilibrium pressure peq = Peq(e, • • • ,?^r^), i-e., it is the 
pressure as computed from the equation of state for given values of e, m , . . . , n„. 
In a general non-equilibrium (dissipative) situation, however, peq is different 
from the isotropic pressure p defined through (15). Denote the difference by 
n = peq — p. Then, the requirement that deviations from local thermodynamical 
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equilibrium are small is equivalent to requiring v^, q^,n^", and 11 to be small 
compared to e, peq, and n,. 

I first outline the phenomenological approach to derive the equations of dis- 
sipative fluid dynamics. For the sake of definiteness, in the remainder of this 
subsection let us consider a system of one particle species only and let us as- 
sume that the total particle number of this species is conserved (implying that 
no annihilation or creation processes take place, i.e., wc do not consider the 
corresponding antiparticles). The particle number current then replaces the net 
charge current. We shall also work in the Eckart frame, where v^^ = 0. Let us 
make an Ansatz for the entropy 4-currcnt S*^. In the limit of vanishing q^^, tt^^ , 
and n, the entropy 4-current should reduce to the one of ideal fluid dynam- 
ics, 5** su^. The only non- vanishing 4- vector which can be formed from the 
available tensors , g^, and ■k'^" is jSq'^, where (3 is an arbitrary coefficient 
(remember n^'^Ui, = 0). Therefore, 



S^ = su^ + pq^' . 



(37) 



With this Ansatz one computes with the help of Ui^d^T^^'^ = and d ■ N = 
h + ne = 0: 



Td- S ={Tp-l)d-q + q-{u + Tdp) + n^^d^u^ + n9>0 
The simplest way to ensure this inequality is to choose 

/3 = 1/T , 

n = ce , 

q^" = kT A^"" {d^lnT - u^) , 



(38) 



(39) 
(40) 
(41) 

(42) 



where (, t], and k are the (positive) bulk viscosity, shear viscosity and ther- 
mal conductivity coefficients. Note that these equations define the dissipative 
corrections as algebraic functions of gradients of the flow velocity u'^ and the 
equilibrium temperature T. With these choices, 



d-S = 



q ■ q 



CT kT^ 



2f]T 



(43) 



which is obviously larger or equal to zero (remember that q ■ q < 0, which can 
be most easily proven from g • u = in the frame where = (1,0)). While this 
ensures the second law of thermodynamics, it was shown [7] that the resulting 
equations of motion are unstable under perturbations and support acausal, i.e., 
superluminous propagation of information. They are therefore not suitable as 
candidates for a relativistic theory of dissipative fluid dynamics. 

A solution to this dilemma was presented by Miiller [8] , and Israel and Stew- 
art [9]. They observed that the Ansatz (37) for the entropy current should not 
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only contain first order terms in the dissipative corrections, but also second order 
terms: 

5" = su'' + /3g'' + Q'' , (44) 

where 

Q^' = ao^q^'+ tt^-^ + u^' {po n'^ + (iiq ■ q + 02 n^^n^x) (45) 

is second order in the dissipative quantities 7T, q^, and tt'"'. Inserting this into 
d ■ S > leads to differential equations for 77, q^, and w^" which involve the 
coefficients C, rj, k, ao, ai, Pn, /3i, /32- It can be shown that, for reasonable values 
of these coefficents, the resulting 14 equations of motion (the 9 equations that 
determine 77, 5'*, and ■k'^" and the 5 conservation equations for N'^, T^") are 
stable and causal. 

In the phenomenological approach, the values of these coefficients are not 
determined. In the second approach, however, based on kinetic theory, they can 
be explicitly computed along with deriving the additional 9 equations of motion 
for 71, qf^, and tt'^^. This will be outlined in the following. 

Let us start by writing the single-particle phase space distribution in local 
equilibrium (20) as 

fo{k,x) = [exp{yo{k,x)} ± 1]"^ , (46) 

where yo{k,x) = [k ■ u{x) — ijl{x)\/T{x). Now assume that the non- equilibrium 
phase space distribution, written in the form 

f{k, x) = ^ [exp{2/(fc, x)} ±1]-' , (47) 

deviates only slightly from the equilibrium distribution function fo{k,x), or in 
other words: 

y{k,x)':^yo{k,x)+ei{x) + k-e2{x) + k^k„s^''{x) , (48) 

where ei{x), st^ix), and £3^ are small compared to T{x), fj,{x). Then one can 
expand f{k,x) around fo{k,x) to first order in these small quantities: 



/(fc,x)~/o(fc,ar) 1 



[y{k,x)-yo{k,x)]] . (49) 



lT^/o(fc,x) 
9 

Note that f{k,x) depends on the 14 variables fi/T — ei, u^/T + £3, and e'^" . 
(£3^" is a symmetric tensor of rank 2, and therefore naively has 10 independent 
components. However, its trace can be; absorbed in a redefinition of the first 
variable fx/T — ei, therefore it actually has only 9 independent components.) 

Inserting f{k,x) into the kinetic theory definition of and T^", (21) 
and (22), (with /o replaced by / and, since we do not consider antiparticles, 
discarding /o), one can establish relations between the 14 unknown macro- 
scopic functions (in the Eckart frame) e, n, u^, 77, g^, tt'"^ and the 14 variables 
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fx/T — ei, uf^/T + ei^, £3^. This uniquely determines the non-equiUbrium single- 
particle phase space distribution f{k,x) in terms of the macroscopic, i.e., fluid- 
dynamical variables. This identification involves one subtlety: as in ideal fluid 

dynamics one still has to know the value of the (equilibrium) pressure to de- 
termine all unknown quantities. The equilibrium pressure Peq iS) however, only 
known as a function of the equilibrium energy density eo and the equilibrium 
particle number density no, but not as function of the actual energy density e 
and particle number density n. Two additional assumptions are required, namely 
that 

e = u^T^^^Uv = eo = u^Tl^^u^ , (50) 
n = u- N = no = u-N() , (51) 

where Tq'^ and Nq are the (kinetic) energy-momentum tensor and particle num- 
ber current computed with the equilibrium phase space distribution fo{k,x). 
Then Pcq{(, n) = Pcq(eoi i^o) a-nd the value of the equilibrium pressure Peq is also 
determined. On close inspection, these additional assumptions do not pose any 
further restriction on the set of 14 unknown functions, but merely serve as defini- 
tions of (equilibrium) temperature T and chemical potential fi corresponding to 
a given energy density e and particle number density n. Another way to say this 
is that the assumptions (50), (51) determine a local equilibrium phase space dis- 
tribution fa{k,x). However, in a non-equilibrium context this distribution has 
no actual dynamical meaning, and one is therefore free to choose it in a way 
which fulfills (50) and (51). 

The next step consists of deriving the equations of motion for the 14 unknown 
functions e, n, u^, U, q^, tt^''. To this end, one takes the first three moments of 
the Boltzmann equation for f{k,x), 

k-df{k,x)=C[f] . (52) 

This results in 

/d^k r H^k 

— k.df{k,x)^d.N = J -^C[/]^0 , (53) 

/d^k f d^k 

— k^k'' fik, x) EE d^T^^^ = y _ fc- C[/] ^ , (54) 

/d'^k f d^k 

— fc^fc^fc^ fik, x) = d^s''-^ = 7 ^"^^^ ^[/] = ■ (55) 

Note that conservation of particle number, energy, and momentum leads to van- 
ishing right-hand sides for eqs. (53) and (54). The structure of the microscopic 
collision term C is such that these requirements are fulfilled (particle number 
and energy-momentum conservation in microscopic collisions between particles) 
[6]. On the other hand, the right-hand side of Eq. (55) does not vanish, since 
there is no corresponding microscopic conservation law. Note that the trace of 
(55) is equivalent to rr? times Eq. (53), such that = 0. Therefore, only 9 out 
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of the set of 10 equations (55) are independent. Together with the 5 equations 
(53) and (54), these 9 equations determine the set of 14 unknown functions of 
dissipative fluid dynamics. The 9 independent equations (55) are equivalent to 
the 9 equations derived from d ■ S > in the phenomenological approach. The 
unknown phenomenological coefficcnts k, tj, ao, ai, /3o, /3i, and /32 can now 
be explicitly identified from suitable projections of X"^. Israel and Stewart have 
shown [9] that the resulting equations fulfill the requirements of hyperbolicity 
and causality. 

This concludes the brief survey of dissipative fluid dynamics. So far, no serious 
attempt has been made to apply rclativistic dissipative fluid dynamics towards 
the description of heavy-ion collisions. First steps were done by Mornas and 
Ornik [10] who investigated the broadening of collisional shock waves through 
dissipative effects in a simple; onc-dimcnsional geometry. Also, Prakash et al. 
generalized the Israel-Stewart theory to a mixture of several particle species 
[11]. 



2.5 Multi-fluid Dynamics 

In multi-fluid dynamics one considers not a single, but several fluids j = 1, . . . , M, 
characterized by the net charge currents N-^j (the net current of conserved charge 
i in fluid j) and energy- momentum tensors Tj^" . There is overall net charge and 
energy-momentum conservation, 

M 



d-Ni=0 , Ni: = ^Ntj , (56) 
i=i 

M 

but not for each fluid separately, 

d ■ Ni, = 5y , a^T/-^ = . (58) 

The right-hand sides deflne the so-called source terms which according to (56), 
(57) obey 

M M 

Y,Sij=0 , ^5;=0 . (59) 

The source terms are parameters of a particular model and have to be specifled 
e.g. from kinetic theory. Let us consider the Boltzmann equation for particles 
from fluid j: 

klra 

The right-hand side involves the collision terms for the microscopic 2-particle 
reactions Im jk (the gain term Cf^) where particles from fluid I and fluid 
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TO {I and TO not necessarily different) collide to produce particles of fluid j and 
k (again, j and k not necessarily different), and jk — > Im (the loss term Cj^) 
where particles from fluid j and k collide to produce particles of fluid / and to. 
Taking the zeroth and flrst moment of this equation yields 

E Sij , (61) 
= . (62) 

This defines the source terms through the microscopic collision rates. 

Results of any specific multi-fluid model will not be discussed here, I in- 
stead refer the reader to the literature on this subject [12]. I close with two 
remarks: (a) a single fluid may consist of several different particle species (for 
instance, tt, K, N, A etc.), as long as it is reasonable to assume that they stay 
in local thermodynamical equilibrium among each other. Then, the only place 
where information enters about these different particle species is the equation 
of state p(e, ni, . . . ,n„). (b) Different fluids may consist of the same particle 
species (with the same equation of state p{e, ui, . . . , n„)). This situation occurs 
for instance in the initial stage of relativistic heavy-ion collisions, where the 
single-particle phase space distributions of target and projectile nucleons, while 
overlapping in space-time, are still well separated in momentum space due to 
the high initial relative velocity between them. This is a situation where there is 
local thermodynamical equilibrium in target and projectile separately, but not 
between them. It therefore is reasonable to treat target and projectile, although 
consisting of the same particle species, as two separate fluids. 



d-Nij = qi j ■dfj{k,x) = QiJ^ 



d^k 



E 



d^k 



k^k''d^f,ik,x) 



kin 



klm 



d^k 



ijk 



nlin 



E 



ijk 



nlm 



3 Numerical Aspects 

In this section, I discuss basic aspects of numerical solution schemes for rela- 
tivistic ideal fluid dynamics. For the sake of simplicity, let us consider the case 
of one conserved charge only. Define 

R = N° =nu° = n-f , (63) 
E = T^^ = {e+ph^-p , (64) 
M^{T0^},^^_^_^ = (e+p)7^v , (65) 

where u^^ = 7(1, v) is the fluid 4-velocity, 7 = (1 — v^)"-*^/^. With these defini- 
tions, the conservation laws (7), (8) take the form 

d- N = dtR + V ■ (i?v) = , (66) 
T^o = dtE + V- [{E + p)v] = , (67) 
T"' = dtM' + V • (MV) + dip = . (68) 
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In this form, the conservation equations can be solved numerically with any 
scheme that also solves the non-relativistic conservation equations. There is, 
however, one fundamental difference between the non-relativistic equations and 
the relativistic ones. In order to solve the latter for R, E, M, the net charge 
density, energy density, and momentum density in the calculational frame, one 
has to know the equation of state p{e, n) and v. The equation of state, however, 
depends on n, e, the not charge density and energy density in the rest frame of 
the fluid. One therefore has to locally transform from the calculational frame 
to the rest frame of the fluid in order to extract n, e, v from R, E, M. In the 
non-relativistic limit, there is no difference between n and i?, or e and E and 
the equation of state can be employed directly in the conservation equations. 
Also, the momentum density of the fluid is related to the fluid velocity by a sim- 
ple expression. The transformation between rest frame and calculational frame 
quantities is described explicitly in the next subsection. 



3.1 Transformation between Calculation Frame and Fluid Rest 
Frame 

In principle, the transformation is explicitly given by equations (63) - (65), i.e., 
by finding the roots of a set of 5 nonlinear equations (the non-linearity enters 
through the equation of state p(e,n)). In numerical applications, however, this 
transformation has to be done several times in each time step and each cell. It 
is therefore advisable to reduce the complexity of the transformation problem. 
This is done as follows [13]. 

First note that M and v are parallel, thus 

M-v = Mw = (e-hp)7^u^ = (e-|-p)(72-l)=£;-e , (69) 

where M = |M|, v= |v|. Therefore, 

e^E-Mv, n = RVl-v^ , (70) 

where the second equation is a simple consequence of (63). With these equations 
e and n can be expressed in terms of R, E, M and v. The 5-dimensional root 
search is therefore reduced to finding the modulus of v for given R, E, and M, 
which is a simple one-dimensional problem. To solve this, use the definition of 
M, 

M = {e+ p)'y^v = {E+ p)v . (71) 
This equation can be rewritten as a fixed point equation for v for given R, E, M: 



E + p{E-Mv,RVl-v'^) 

The fixed point yields the modulus of the fiuid velocity, from which one can 
reconstruct v = t; M/M, and find e and n via (70). The equation of state p(e, n) 
then yields the final unknown variable, the pressure p. 
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3.2 Operator Splitting Method 

In general, to model a heavy-ion collision with ideal fluid dynamics requires 
to solve the 5 conservation equations in three space dimensions. Since this is in 
general a formidable numerical task, one usually resorts to the so-called operator 
splitting m,ethod, i.e., the full 3-dimensional solution is constructed by solving 
sequentially three one- dimensional problems. More explicitly, all conservation 
equations are of the type 

dtU+ diFi{U)=0 , (73) 

i=x,y,z 

U being R, E, or M*. Such an equation is numerically solved on a space- time 
grid, and time and space derivatives are replaced by finite differences: 

Ui;f = Urjk-AtG[Urjk\ , (74) 

where i, j, k are cell indices (the cell number in x, y, and z direction) and n 

denotes the time step. At is the time step width. G U-j^. is a suitable finite 

difference form of the 3-divcrgcncc in (73). 

It can be shown that in the continuum limit instead of solving (74) it is 
equivalent to solve the following set of predictor-corrector equations 

^ijk "^"^ = U^jk - Gx [Uiji^ , 

(75) 



uSl^^' = U^^'-AtG, 



7-r(l)"+l 
'-'ijk 

rr(2) n+1 



(76) 



and that the solution converges towards the solution of (73). Here, the Gi[U], 
i ~ x,y,z, arc finite difference forms of the partial derivatives di Fi{U) (no 
summation over i) in x, y, or z direction, t/j^^"^^ is the first prediction for 
the full solution U^jl^- It is generated by solving a finite difference form of the 
one-dimensional equation 

dtU + diFi{U) = , (77) 

where i = x. Subsequently, the first prediction uj^jl """"^ is used to solve a finite 

difference form of (77), where now i = y,to obtain the second prediction uj-jl """"^ 

for the full solution. {U^l"'^^ has been corrected to u\Jl'^~^^ .) Finally, the full 

solution U^p.^ is obtained by using U-jl"^^ to solve a finite difference form of 

(77) with i = z. {U^jl''^^ has been corrected to J7^+\) 

In other words, the solution to the partial differential equation (73) in three 
space dimensions is obtained by solving a sequence of partial differential equa- 
tions (77) in one space dimension. The 3-divergence operator in (73) was split 
into a sequence of three partial derivative operators. Physically speaking, in a 
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given time step one first propagates the fields in x direction, then in y direc- 
tion, and then in z direction. In actual numerical applications it is advisable to 
permutate the order xyz to minimize systematical errors. 

The advantage of the operator splitting method is that there exists a vari- 
ety of numerical algorithms which solve evolution equations of the type (77) in 
one space dimension (see, for instance, [14] and refs. therein). One of them is 
discussed in the following subsection. 

3.3 The Relativistic Harten Lax van Leer Einfeldt Algorithm 

The relativistic Harten-Lax-van Leer-Einfeldt (HLLE) algorithm [14, 15] solves 
equations of the type 

dtU + d,F{U) = , (78) 

i.e., propagation of a field U in one space dimension. For ideal relativistic fluid 
dynamics, U = R, E, ov M and F{U) = Rv, {E + p)v, or Mv + P- (For one- 
dimensional propagation, it is sufficient to consider only the components of the 
momentum density M and the fluid velocity v in the direction of propagation. 
They are here denoted by M and v, respectively.) 

The idea behind the relativistic HLLE scheme is the following. Consider the 
initial distribution of the density [/ on a numerical grid. U is assumed to be 
constant inside each cell, but different from cell to cell, i.e., the initial distribu- 
tion consists of a sequence of constant flow fields inside the cells separated by 
discontinuities at the cell boundaries, cf. Fig. 1. 



Ui.2 




i-2 i-1 i i+1 i+2 i+3 

< > 

Ax 

Fig. 1. The initial distribution of the density U on the numerical grid. 

In the further time evolution these discontinuities will decay, resulting in the 
transport of U across the grid. The decay of a discontinuity between two regions 
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of constant flow is, however, a well-known problem in fluid dynamics, the so- 
called Riemann problem. For simple equations of state it is even analytically 
solvable. Consider the discontinuity to be located at x = 0. Denote the density 

in the region of constant flow to the left of the discontinuity by Ui, and that to 
the right by 1/^. The initial condition at time t — then reads 

cf. Fig. 2 (a). For the sake of definiteness, consider Ui > U^. For t > 0, the solu- 
tion looks qualitatively as in Fig. 2 (b). There is a rarefaction fan propagating 
into the region of higher density (in this case to the left), and a shock front into 
the region of lower density (in this case to the right). Between fan and shock 
wave there are two regions of constant flow separated by a contact discontinuity 
(a discontinuity where the pressure is equal on both sides). It is evident that a 
numerical algorithm can be constructed which solves the fluid dynamical equa- 
tions simply by solving a sequence of Riemann problems for the discontinuities 
at all cell boundaries in a given time step. Such algorithms are called Godunov 
algorithms [16]. 

The relativistic HLLE is a so-called Godunov-type algorithm [16], i.e., it does 

not employ the full solution of the Riemann problem but approximates it by a 
region of constant flow between Ui and Ur, cf. Fig. 2 (c): 

X < bit 

U{x,t) = { Uir , bit<x <brt . (80) 

X>bi:t 

Here, 6i < and &r > are the so-called signal velocities. They character- 
ize the velocities with which information about the decay of the discontinuity 
travels to the left and right into the regions of constant flow. The value Uir 
in the region of constant flow between Ui and Ui- is determined in accordance 
with the conservation laws. To this end, integrate (78) over a flxed interval 
[s^minj 2;niax]) Xjnin < bit, Xmax > b^t. One obtains: 

b,Ur-biUi + F{Ui)-F{U,) 

The value of the flux F(Uii) corresponding to the density Uir is determined by 
integrating (78) over the fixed interval [0, aimax] or [a;min,0]: 

= b.FiUi)-biFiU.) + bib.iU.-Ui) 

br - h 

Upon discretization, the differential operator dx F{U) in the evolution equation 
for the density Ui in cell i assumes the form [F(C/i+i/2) — F{Ui-i/2)\/ Ax where 
Ax is the cell size (grid spacing) and Ui±i/2 are values of the density at the 
position of the right and left boundary of cell i. These values are taken after 
the decay of the respective discontinuities at the cell boundaries, i.e., they are 
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(a) 



t=0 



^ X 







u, 



rarefaction fan 



(b) 




contact discontinuity 
shock front 



^ X 



Ui 




— ' ■ — > X 

Fig. 2. (a) The initial condition of the Riemann problem at t — 0. (b) The solution 
of the Riemann problem at t > 0. (c) The approximate solution of a Godunov-type 
algorithm. 



the corresponding values Uir given by (81) and the respective -^'(£^^±1/2) are 
the corresponding vahios F{Uii) given by (82). This yields the following explicit 
expressions for the relativistic HLLE scheme 



n+l 



Ax 



F 



br F - 61 F ([/f+o + br 61 {ur+1 - ur) 

br - h 



(83) 
(84) 



A reasonable estimate for the signal velocities is to take them as the relativistic 
addition (subtraction) of flow velocities and sound velocities in the respective 
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cells adjacent to the cell boundary: 

6. = max|o, ^j-;+'^+^ 1 , (85) 

As described above, this scheme is accurate to first order in time. A scheme 
which is accurate to second order can be obtained using half-step updated values 
F {Ui±x/2^^ more details see [17]. 



4 One-dimensional Solutions 



In this section I discuss solutions of ideal relativistic fluid dynamics in one 
space dimension. 1 first introduce the notion of characteristic curves. Then, 1 
discuss possible one-dimensional wave patterns for thermodynamically normal 
and anomalous media. Choosing a representative equation of state which fea- 
tures both thermodynamically normal and anomalous regions I then discuss the 
expansion of semi-infinite matter into the vacuum. The emerging wave patterns 
will help us to miderstand the possible solutions of the Landau model, which 
was historically the first fiuid-dynamical model for relativistic heavy-ion colli- 
sions. Finally, also the Bjorken model for ultrarelativistic heavy-ion collisions is 
discussed. 



4.1 One-dimensional Wave Patterns 

For flow in one spatial dimension (say, in x direction) the two conservation 
equations for energy and for momentum read: 

dt T°° -h = , dt -h T^^ = . (87) 

A suitable linear combination of these equations leads to the equivalent set of 
equations 

'dt + ^^d^n±^i) , (88) 

l±VCs J 

where Cg = dp/de\s/n is the velocity of sound squared (s/n is the specific entropy) 

TZ±^y-yo± [ de' ^''f (89) 
Je„ e'+P(e') 

are the so-called Riemann invariants, y = Artanhw is the fluid rapidity. Equation 
(88) has the obvious interpretation that the Riemann invariants TZ± are constant 
along world lines x± {t) defined by 

^^w^ = ^^. (90) 
dt l±vcs ^ ' 
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These world lines are the so-called characteristic curves or characteristics C± {x,t). 
It is also obvious that these curves are the world lines of sonic perturbations or 
sound waves on top of the fluid- dynamical wave pattern. C+{x,t) characterizes 
sound waves moving to the right (in positive x direction) while C_ (x, t) charac- 
terizes those moving to the left (in negative x direction). For the simple example 
of constant flow, the characteristic curves are shown in Fig. 3. 




Fig. 3. The characteristic curves for a constant flow pattern. 



Let us now consider a so-called simple rarefaction wave moving to the right, 
cf. Fig. 4. (For the definition of a simple wave, see [18], for our purposes it is suf- 
ficient to remark that in one spatial dimension a simple wave is the only possible 
wave that can connect two regions of constant flow. A rarefaction wave denotes a 
wave where the energy density decreases in the direction of propagation.) Then, 
one can prove that TZ+ = const, everywhere (for the proof, see [18]; analogously, 
for simple waves moving to the left, TZ- = const.). 




Fig. 4. A continuous simple wave between two regions of constant flow, moving to the 
right. 



It is therefore suSicient to consider the equation for the TZ- invariants, or 
the C- characteristics, rcspcctivcily. Let us consider how the slope of the C_ 
characteristics changes with x at constant t: 



dx 



w_ 



(91) 
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Prom TZ+ = const, everywhere one infers 



while 



Therefore, 



where 



V = —(1 — V ) e , 

^ ^ e + p 



(92) 
(93) 

(94) 
(95) 



Equation (94) is an important quahtative result: Since the first factor is always 
positive {w- as well as Cg are causal), and since the energy density decreases 
with X for the rarefaction wave considered here, e' < 0, the sign of w'_ is solely 
determined by the sign of S. The quantity E, however, is solely determined by 
the equation of state of matter under consideration, i.e., its sign (and absolute 
value) is an intrinsic property of the fluid. Matter with Z" > is called thermo- 
dynamically normal, while matter with U < is therm,odynamically anomalous. 
More specifically, ii E > 0, then w'_ > 0, and ii E < 0, then w'_ < 0. A positive 
w'_, however, means that the C_ characteristics "fan out" in the x — t plane, 
while a negative w'_ indicates that they converge and ultimately intersect at one 
point, cf. Fig. 5. 



t t 




(a) E>0 >0 (b) E<0 \^'<0 

Fig. 5. For a simple wave moving to the right and (a) 17 > the C- characteristics fan 
out, while for (b) X" < they converge and intersect. 
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Intersecting characteristics, however, signal the formation of shock waves. 
Physically speaking, picture a sonic perturbation (travelling along a character- 
istic) emitted at a point Xi. This perturbation will eventually overtake a per- 
turbation emitted at X2 > xi (namely when the corresponding characteristics 
intersect). Thus, the two small perturbations add up to form a larger one. Imag- 
ine this happening for other perturbations (emitted at different points) as well. 
Eventually, a finite discontinuity (shock wave) is formed from the siiperposition 
of a large number of infinitesimal sonic perturbations. Shock waves are discon- 
tinuous solutions of ideal fluid dynamics and will be discussed in more detail in 
the following subsection. 

I conclude this subsection by collecting the above arguments in the following 
classification scheme of one-dimensional wave patterns. Continuous rarefaction 
waves arc stable in therm,odynamieally norm,al matter while they are unstable 
in anomalous matter. On the other hand, rarefaction shock waves are stable in 
thermodynamically anomalous matter while they are unstable in thermodynam- 
ically normal m,a,tter. If we perform an analogous consideration for a continuous 
compression wave we are led to the conclusion that such waves are unstable in 
normal and stable in anomalous matter, while compression shock waves are sta- 
ble in normal and unstable in anomalous matter. These results are summarized 
in Table 4.1. A "+" sign means "stable" while a "— " sign indicates "unstable". 



Table 1. Classification scheme for the stability of one-dimensional wave patterns. 



Wave 



r > 



r < 



Continuous rarefaction 

Rarefaction shock 
Continuous compression 
Compression shock 



+ 
+ 



Most matter is thermodynamically normal. In the presence of phase transi- 
tions, however, an equation of state can feature regions where matter is ther- 
modynamically anomalous. As will be seen in Subsections 4.4 and 4.5, this will 
strongly influence the time evolution of the system in a qualitative and quanti- 
tative way. 

4.2 Shock discontinuities 

Shock waves represent discontinuous solutions of ideal fluid dynamics. While the 
partial derivatives of N^ and T^^" appearing in the conservation equations are ill- 
defined at the location of such discontinuities, there is still a simple way solve the 
problem of charge and energy-momentum transport across a shock discontinuity. 
To this end, let us consider the case of one conserved charge only, and study such 



Fluid Dynamics for Relativistic Nuclear Collisions 



23 



a discontinuity in its rest frame. Matter enters the discontinuity with velocity vq 
in a thermodynamic state characterized by the net charge density no, the energy 
density cq, and the pressure po (which is of course determined by eo and no 
through the equation of state). The task is to determine the velocity v and the 
thermodynamic state of matter (n, e, and p) emerging from the shock. Imagine 
a small volume V which encloses the discontinuity, cf. Fig. 6. 



£() Po "() 



£ p n 



Fig. 6. A shock discontinuity in its rest frame. 



Let us now integrate the conservation equations (7) (for a single conserved 
charge) and (8) for one-dimensional flow over V: 



Jv Jv 



a / d^xTOO 



d^xa^T^" = , 



dt I d^xT^°+ / d^xa^T^^ = 
Jv Jv 



(96) 
(97) 
(98) 



In a steady-state situation (a stable, propagating shock discontinuity) the total 
amount of charge, energy and momentum inside V cannot change with time, 
therefore, the first terms in these equations vanish. The other terms are inte- 
grated by parts to yield the set of equations 



nolot'o 

{e-o + Po)llvl + Po 



(e+p)j'^v , 
(e +p)^'^v'^+p 



(99) 
(100) 
(101) 



These are the conservation equations for net charge and energy-momentum 
across a shock discontinuity. They are no longer partial differential equations, 
but purely algebraic. For a given initial state no, cq, Po, and velocity vq, they de- 
termine the final state n, e, p, and the velocity v of compressed matter emerging 
from the shock, if the equation of state p{e, n) is known. 
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One can eliminate the velocities from the set of equations (99) - (101) to 
obtain the so-called Taub equation [19] 

(e + p)X - (eo + Po)Xo = {p - po){X + Xo) , (102) 

where X = {e +p)/'n? is the so-called generalized volume. Once p(c, n) is fixed, 
the solution of the Taub equation defines the so-called Taub adiabat p{X), cf. 
Fig. 7. For a given initial state {po,Xo) (the so-called center of the adiabat) it 
represents all final states {p, X) for matter emerging from the shock, which arc in 
agreement with net charge and energy-momentum conservation. The actual final 
state is then selected by specifying vq . This determines all variables uniquely in 
the rest frame of the shock. The remaining unknown is, however, the velocity 
of the shock in an arbitrary calculational frame. For compressional shock waves, 
such as occur in the initial stage of heavy-ion collisions (cf. [20] for a detailed 
discussion), this shock velocity can be uniquely determined from the geometry 
of the collision. For rarefaction shock waves this is not possible, and thus in 
principle there is a whole region of final states on the Taub adiabat, which are 
in agreement with energy-momentum and net charge conservation. It turns out, 
however, that the stationary situation is always given by a rarefaction shock 
where matter emerges at the so-called Chapman- J ouguet point, indicated by 
"C J" in Fig. 7 (b) [5] . This point is defined as the point where a chord between 
the center (po, -'^o) and a final state on the adiabat is tangential to the adiabat. 
This then uniquely fixes the state of matter emerging from the shock, as well as 
the velocity of the shock in the calculational frame. Note that it is also possible 
to define a Taub adiabat in the case that there is no conserved charge, see [17, 21] 
for details. 

To conclude this subsection, let us consider what happens to the entropy 
flux across a shock discontinuity. Integrate the conservation equation (35) for 
the entropy current over the volume V which encloses the shock front in its rest 
frame, 

dt [ d^xs7-F / d^xd^s'yv = , (103) 
Jv Jv 

and perform an integration by parts in the second term. This yields: 

s^v = SQ-^QVQ + -^dtS , (104) 

where Ai_ is the transverse area of the shock front and S = /yd^a:;s7 is the 
total entropy inside the volume V . The second law of thermodynamics tells us 
that the entropy cannot decrease, dtS >Q. Consequently, 

s^v> sojovo ■ (105) 
Dividing both sides by (99) one concludes 



(106) 
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P P 




Fig. 7. (a) The Taxib adiabat for a comprcssional shock wave. {pq,Xo) is the center 
of the adiabat, {p,X) is one final state on the adiabat which is selected by a choice 
of vo- (b) The Taub adiabat for a rarefaction shock wave. "CJ" denotes the Chap- 
man- Jouguet point. 

i.e., the specific entropy increases across a shock front. This result is remarkable, 
since we know that the entropy current is conserved in ideal fluid dynamics, 
Eq. (35). However, this equation holds strictly only for continuous (differen- 
tiable) solutions. Shock discontinuities do not belong to this class, and therefore 
can produce entropy. Physically speaking, microscopic non-equilibrium processes 
take place inside a shock front which lead to this increase of entropy. 

One could object that this conclusion is not stringent in the sense that (106) 
also allows for the case where s/n = sq/uq, i.e., where the entropy does not in- 
crease across the shock front. However, by explicitly solving the shock equations 
(99) - (101) with a given equation of state one finds that this case occurs only 
for infinitesimal shock discontinuities (which then degenerate into sonic pertur- 
bations, which in turn preserve entropy). For any finite discontinuity one finds 
s/n > So/no. 

The Chapman-Jouguet point (cf. Fig. 7) is actually special in this respect: 
it corresponds to that state of matter emerging from the shock, where entropy 
production is maximized [5]. It is amusing to note that in selecting this state 
as the final state of matter emerging from a rarefaction shock wave (cf. discus- 
sion above), fluid dynamics not only automatically respects the second law of 
thermodynamics, but even exploits it to the maximum extent. 

4.3 Equation of State and Expansion into Vacuum 

In this subsection I discuss possible wave patterns for the one-dimensional ex- 
pansion of semi-infinite matter into the vacuum. To be specific, let us first choose 
an equation of state which bears relevance to relativistic heavy-ion physics. At 
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zero net baryon number, QCD lattice data [4] suggest the following Ansatz for 
the entropy density as function of temperature: 

^ ^^^3l-tanh[(T-Te)/^T] ^ ,^^3 1 + tanh[(T- T.)/^T] ^^^^^ 

where cq/ch is the ratio of degrees of freedom in the quark-gluon phase and 
the hadronic phase, Tc ~ 160 MeV is the (phase) transition temperature, and 
AT is the width of the transition. Present lattice data are not yet sufficiently 
precise to decide whether the transition is first (corresponding to AT = 0) or 
higher order, or just a smooth cross-over transition, but they restrict AT to 
be within the range < AT <, 0.1 Tc. Note that for AT = the equation of 
state becomes that of the well-known MIT bag model [22] with a bag constant 
B = (cq/ch — l)Pc) where Pc is the pressure at the phase transition temperature 
T 

To cover the possible range of AT, we shall consider the limiting values 
AT = and AT = 0.1 T,, in the following. Both cases will be compared to 
results for an equation of state where there is no transition to the quark-gluon 
phase, i.e., where 

s{T) = sh(T) = chT^ . (108) 

Once s{T) is known one can compute other thermodynamic variables from fun- 
damental thermodynamic relations, for instance: 



nT 

P 



[ dT's{T') , e = Ts-p . (109) 
Jo 



The three equations of state considered here are explicitly shown in Fig. 8. The 
ratio of degrees of freedom cq/ch was chosen to be 37/3, corresponding to an 
ultrarelativistic gas of u and d quarks and gluons in the quark-gluon phase and 
a massless pion gas in the hadronic phase. 

Figs. 8 (a,b) show the entropy density divided by T^ and the energy density 
divided by as fiuictions of T. This representation of the equation of state is 
commonly used by the lattice QCD community, On the other hand, fluid dynam- 
ics requires the pressure as a function of energy density, p(e) , which is shown in 
Fig. 8 (c). The collective evolution of the fluid is, however, controlled by pressure 
gradients. Figure 8 (d) shows the velocity of sound squared Cg = dp/de (if there 
are no conserved charges). This quantity determines the pressure gradient dp for 
a given gradient in energy density de, i.e., it characterizes the capability of the 
fluid to perform mechanical work, or in other words, it characterizes the expan- 
sion tendency. Thus, for the equation of state with a first order phase transition, 
AT = 0, in the mixed phase of quark-gluon and hadronic matter, fn < c < cq, 
the system does not perform mechanical work and therefore has no tendency to 
expand. As will be seen in the following this will have profound influence on the 
time evolution of the system. 

For the equation of state with a smooth cross-over transition, AT = 0.1 Tc, 
the expansion tendency is not zero, but still greatly reduced in the transition 
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Fig. 8. (a) The entropy density divided by as a function of T. (b) The energy density 
divided by as a function of T. (c) The pressure as a function of energy density, (d) 

The velocity of sound squared as a function of energy density, cq/ch = 37/3. Units of 
energy are Tc, units of energy density are T^Sc, where Sc = (cq + ch) T'c /2. Solid line: 
AT = 0, dotted line: AT = 0.1 Tc, dashed Une: ideal hadron gas. 



region as compared to the ideal gas equation of state without any transition 
(Cg = 1/3 = const, for all values of e). The transition region en ^ e ^ eg is 
referred to as the "soft region" of the equation of state [23] . For an equation of 
state with a first order transition, the point e = eg is called the "softest point" 
of the equation of state [24] . (This notion comes from considering the function 
p(e)/e which has a minimum at eq.) 

Another quantity of interest is E, which determines whether matter is ther- 
modynamically normal or anomalous. Figure 9 shows this quantity (times Ts) 
as computed from (95) for the three equations of state studied here. For AT = 0, 
matter becomes anomalous in the mixed phase, the other two equations of state 
are thermodynamically normal everywhere. (Strictly speaking, 17 — only van- 
ishes in the mixed phase, but does not become negative. This is, however, suffi- 
cient for the formation of stable rarefaction shock waves.) 

Let us now consider the onc-dimcnsional expansion of semi-infinite matter 
into the vacuum. Figure 10 shows temperature profiles for (a) the expansion of 
an ideal gas and (b,c) for the expansion with the AT = equation of state. In 
(b) the initial energy density of semi-infinite matter is chosen to be well above 
eg, the phase boundary between the quark-gluon and the mixed phase, in (c) the 
initial energy density is just below eg. The dotted line in (a) indicates the initial 
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Fig. 9. The quantity U (times Ts) as a function of e for AT = (solid line), 
AT = 0.1 Tc (dotted line), and the ideal hadron gas equation of state (dashed line). 



temperature profile for all cases. The initial profile indicates a discontinuity 

at .T = which separates two regions of constant flow, the semi-infinite slab 
of matter at rest to the left {x < 0), and the vacuum to the right {x > 0). 
This initial condition is in fact a special case of the Riemann problem discussed 
in Subsection 3.3. From general arguments (see above) the solution at f > 
can only be a simple wave, connecting these two regions of constant flow. For 
the ideal hadron gas which is thermodynamically normal matter, we have seen 
above that this simple wave must be a continuous rarefaction wave, in this case 
moving to the right. As mentioned above, for such a wave the Riemann invariant 
TZ+ = const, everywhere, cf. (89), from which we deduce the relationship between 
the fluid rapidity y and the energy density e on the rarefaction wave: 

2/(e) = -T^ln- , (110) 

where we have used the fact that for the ideal hadron gas equation of state 
p = CgC and that the initial fluid rapidity of the semi-infinite slab is zero, yo = 0. 
The fluid velocity on the rarefaction wave is then given by v^e) = tanhy(e). 
Finally, the position at which one finds a given energy density e at time t can 
be deduced by integrating (90) for the non-trivial C_ characteristics: 

where we have used the fact that the initial position of the simple wave is at 
a; = and that Cg = const, for the ideal hadron gas equation of state (we have 
assumed that the hadron gas consists of massless, i.e., ultrarelativistic pions, for 
which Cg = 1/3). Equation (111) tells us that the rarefaction wave moves with 
sound velocity into the semi-infinite slab of matter (to the left), xa = —Cst, and 
with the velocity of light into the vacuum (to the right), xb = t. 
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Fig. 10. Temperature profiles for the expansion of semi-infinite matter into vacuum, 
(a) Ideal hadron gas equation of state, the dotted line indicates the initial state, the 
temperature is normalized to the initial temperature To. (b.c) Equation of state with 
AT = 0, in (b) the initial energy density is well above eg, in (c) it is just below eg. 
The temperature in (b,c) is normalized to the critical temperature Tc. 



The expansion in the case of a first order phase transition, AT = 0, proceeds 
similarly, with the exception that in the region of energy densities corresponding 
to the mixed phase, matter is thermodynamically anomalous, cf. Fig. 9, such 
that from Table 4.1 we conclude that the stable wave pattern is not a continuous 
rarefaction wave, but a rarefaction shock wave. Thus, as long as matter is in the 
(thermodynamically normal) quark-gluon phase, the expansion will proceed as a 
continuous rarefaction wave as in Fig. 10 (a), but upon entering the mixed phase 
(energy density eg, temperature Tc) a rarefaction shock wave will form. The state 
of matter emerging from this shock wave is determined from the shock equations 
as described in the previous subsection, i.e., it corresponds to the Chapman- 
Jouguet point on the Taub adiabat with center located at the phase boundary 
between quark-gluon and mixed phase (for more details, see [17]). Then, also 
the velocity of the shock Vsh in the calculational frame is determined. In general 
Vgh and the velocity of matter at the base of the continuous rarefaction wave 
are not equal. This leads to the formation of a plateau of constant fiow between 
xb and xc- The state of matter at the Chapman- Jouguet point corresponds 
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to thermodynamically normal hadronic matter, so that the further expansion 
has to proceed as a continuous rarefaction wave. The emerging wave pattern is 
shown in Fig. 10 (b). 

The only difference between Fig. 10 (b) and (c) is that the initial energy 
density in (c) is already below eg, i.e., in the region corresponding to mixed 
phase. Therefore, the expansion starts out with a rarefaction shock wave, from 
which matter emerges at the Chapman Jouguct point of the respective Taub 
adiabat with center corresponding to the initial state of matter. (Note that this 
Taub adiabat differs from the one in (b), since their centers are different.) Further 
expansion proceeds as a continuous rarefaction wave in hadronic matter. 

This completes the discussion of the expansion of semi-infinite matter into 
vacuum and prepares us to understand the Landau model which is subject of 
the next subsection. 

4.4 The Landau Model 

The Landau model is historically the first case where fluid dynamics was applied 
to describe at that time hadron-hadron collisions [25]. Its main focus of ap- 
plication nowadays is, of course, nucleus- nucleus collisions. The main ideas are 
summarized in Fig. 1 1 . Imagine two nuclei colliding at ultrarelativistic velocities 
in their center of mass. The nuclei are Lorentz-contracted to a "pancake-like" 
shape. In the moment of impact, nuclear matter becomes highly excited (the 
detailed microscopic processes which happen during this stage are of no concern 
for the following). In the limit that the velocities of the nuclei v 1, there will 
be no baryon stopping (due to the limited stopping power of nuclear matter), 
i.e., the baryon charges will pass through each other unscathed, leaving highly 
excited, net baryon-free matter in their wake. Due to Lorentz contraction, the 
initial extension 2 L in 2: direction of this slab of highly excited matter is much 
smaller than the transverse size of the slab, such that the expansion will pro- 
ceed mainly in the longitudinal direction and is thus essentially one-dimensional. 
The Landau model assumes that the slab has no initial collective velocity and 
that rapid thermalization takes place which is completed at t = 0, It is also 
assumed that the equation of state has the simple ultrarelativistic form p = c^e, 
Cg = const., i.e., that matter is thermodynamically normal for all e. (The orig- 
inal idea of Landau actually was that the baryons are immediately stopped in 
the collision through compressional shock waves. Data from heavy-ion experi- 
ments at BNL-AGS and CERN-SPS prove that this picture is unrealistic, due 
to the aforementioned finite stopping power of nuclear matter. However, since 
the collision is ultrarelativistic, the thermal energy in the highly excited slab 
is much larger than the chemical energy associated with the conservation of 
baryon charge. Therefore, to good approximation, /xb = riB = 0, and the further 
evolution of the slab will be identical to what is discussed here.) 

For t > 0, the slab starts to expand. As in the expansion of semi-infinite 
matter discussed in the previous subsection, rarefaction waves will form. For 
thermodynamically normal matter, these are continuous (Riemann) rarefaction 
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Fig. 11. The Landau model for nuclear collisions. See text for details. 



waves which travel into the slab with sound velocity. Therefore, they will meet 
at the center of the slab (here chosen to be the origin z = 0) at a time i/cg. For 
times t > L/cs, these waves overlap and the solution becomes more complicated. 
In a region near the hght cone, the solution will remain a Riemann rarefaction 
wave, therefore we term this region the Riemann region. In the center where 
the Riemann rarefaction waves overlap, however, the solution is no longer a 
simple wave (indeed, only two regions of constant flow have to be connected 
by a simple wave [18], for two simple waves no such theorem exists). For Cg = 
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const, the solution can still be given in closed analytic form [25], although the 
derivation is rather complicated. However, since two of our equations of state do 
not have constant velocities of sound, we have to resort to numerical solution 
methods, such as the relativistic HLLE discussed above. In principle, numerical 
algorithms can deal with arbitrary (physically reasonable) equations of state, 
and are therefore well able to handle this problem (although one should test 
them thoroughly for test cases where analytical solutions are known [17, 20]). 




Fig. 12. Expansion in the Landau model for AT = (a,d), AT = 0.1 Tc (b,e), and the 
ideal gas equation of state (c,f). (a-c) show temperature profiles for different times, 
(d-f) show the corresponding isotherms in the t — z plane (numbers are temperatures 
in units of Tc). The initial energy density is eo = 1.875 TcSc in all cases. 



In Fig. 12 numerical solutions for the Landau model are presented for the 
three different equations of state of Fig. 8. The initial energy density is eo = 
1.875 TcSc which is slightly larger than eg. In Figs. 12 (a-c) temperature profiles 
are shown for different times t and for the 2; > half plane (the solution in 
the other half plane is the respective mirror image). For AT = 0, Fig. 12 (a), 
one clearly observes the rarefaction shock wave which, for this initial energy 
density is almost stationary. Hadronic matter is expelled from the shock until 
the energy in the interior of the slab decreases below en and the shock vanishes. 
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For AT = 0.1 Tc, Fig. 12 (b), no shock is formed, although the variation of 
the velocity of sound in the mixed phase, Fig. 8 (d), leads to shapes for the 
continuous rarefaction waves which differ strongly from those for a constant 
velocity of sound, Fig. 12 (c). Note the kink in the temperature profiles in the 
latter case which indicate the position where the Landau solution matches to the 
Riemann rarefaction wave. Note also the difference in the initial temperatures 
for the three cases although the initial energy density is the same. This is a 
consequence of the different number of degrees of freedom for the three equations 
of state at high energy densities. 

In Figs. 12 (d f) corresponding isotherms are shown in the t — z plane. The 
most pronounced feature is that due to the small propagation velocity of the 
rarefaction wave, the system stays hot for a much longer time span for the 
AT = equation of state, Fig. 12 (d), than for the ideal gas. Fig. 12 (f). This 
is in agreement with the general argument presented earlier that the softening 
of the equation of state in the mixed phase region leads to a reduced expansion 
tendency and thus to a "stalled" expansion of the system. The softening of the 
equation of state is also the reason why the expansion for the AT = 0.1 T^, 
equation of state. Fig. 12 (e), is delayed in comparison to the ideal gas case, 
although no rarefaction waves are formed. For a quantitative analysis of the 
delayed expansion in the Landau model see [23] . 

4.5 The Bjorken Model 

One of the main assumptions of Landau's model is that the initial collective 
velocity of the slab of excited matter vanishes. However, this cannot be quite 
true on account of the following argument. In the limit u — * 1, the size of the 
nuclei in longitudinal direction goes to zero, and there is no scale in the problem 
at all. In this case, the collective velocity of matter in the slab has to be of the 
scaling form v = z/t. The consequences of this special form for the longitudinal 
fluid velocity were first investigated in [26, 27], again with respect to possible 
applications in hadron-hadron collisions. Bjorken [28] was the first to discuss it 
in the framework of nuclear collisions, 

The main ideas of the so-called Bjorken model are summarized in Fig. 13. 
As in the Landau model, two ultrarelativistic, Lorentz-contracted nuclei collide 
at ^; = and t = (the moment of complete overlap) in the center of mass 
frame of the collision. Due to the limited amount of nuclear stopping power, the 
baryon charges keep on moving along the light cone, while microscopic collision 
processes (the nature of which is of no concern for the following) lead to the 
formation of a region of highly excited, net charge free matter in the wake of 
the nuclei. In contrast to the Landau model, however, the collective velocity in 
this region is of the scaling form v = z/t. The region of highly excited matter 
is supposed to rapidly equilibrate locally within a time span tq (which is of the 
order of a fm or less), and the further evolution of the system can be described in 
terms of ideal fiuid dynamics. One important point is that, due to the absence of 
a scale, physics has to be the same for matter at different longitudinal coordinate 
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z if compared at the same proper time t = t\/l — v'^ = \Jt'^ — . (Such curves of 
constant r describe hyperbola in space-time.) Thus, the initial thermodynamic 
state of all fluid elements is the same at the same proper time tq. 




If the longitudinal velocity profile is enforced by the scaling argument, the 
fluid-dynamical solution simplifies in fact considerably. To see this, change the 
variables z in the conservation laws for one-dimensional longitudinal motion 
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in the absence of conserved charges, 
9tT°° + a^T^° = 



dt T"^ + d, T^^ = 



(112) 



to the variables r = ^t'^ — z^, which is the proper time of a fluid element, and 
r\ = Artanhi; = Artanh[2/f], which is the rapidity of a fluid element. Then, the 
coupled system of partial differential equations (112) decouples: 



de 



+ 



e + p 



dp 
dr] 



. 



(113) 
(114) 



The second equation (114) has the interesting consequence that there is no pres- 
sure gradient between adjacent fluid elements (the one at t] and the one at 
rj + drj). At first glance this would seem to indicate that there is no expansion 
of the fluid at all. This, however, is not true, since the fluid velocity is certainly 
finite, V = z/t. The answer is that the new coordinates (r, rj) already take the 
scaling expansion into account: a fluid element at -q with a width Arj in fact 
"grows" in longitudinal direction by an amount dz = AtArj during the time 
span dt . 

Another consequence of (114) is derived from the Gibbs-Duhem relation: 



dp 
dr] 



dT 
dr] 



Ed]ii 



(115) 



This equation means that for vanishing conserved charges rij = 0, i = 1, . . . , n, 

the temperature has to be constant along curves of constant r, i.e., along the 
space-time hyperbola shown in Fig. 13 (77 varies along these curves). In the 
general case of non-zero net charges, however, only the particular combination 
of charge densities, entropy density, and derivatives of T and the ]ii appearing 
in (115) has to vanish along curves of constant r. Equation (114) represents the 
principle of "boost invariance" commonly associated with the Bjorken model: at 
constant t the pressure is independent of the longitudinal rapidity, i.e., it is the 
same in fluid elements with different 77, or in other words, it does not change 
if one performs a longitudinal boost to a different reference frame. This is a 
consequence of the scaling form for the longitudinal velocity. 

Equation (113) also has an interesting consequence. With the first law of 
thermodynamics, one derives as usual the conservation of the entropy current 
which now takes the form 



ds_ 
d^ 



+ - =0 

T 



(116) 



which can be immediately integrated to give 



ST = sqTo — const. 



(117) 
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at constant 77. The constant may in principle differ for different 77, but since 
the initial thermodynamic state along tq was the same for all i]. that constant 
will also be the same for all rj at other r > tq. Equation (117) is interesting 
because it tells us that the entropy density decreases inversely proportional to r 
independent of the equation of state of the fluid. The time evolution for energy 
density, pressure, or temperature might depend on the equation of state, but not 
the one for the entropy density. 




Fig. 14. Proper time evolution for (a) energy density, (b) entropy density, (c) pressure, 
and (d) temperature in the Bjorkcii model for nuclear collisions (k)ngitndinal expansion 
only). Solid line: AT = 0, dotted line: AT = 0.1 Tc, dashed line: ideal gas equation of 
state. The initial energy density is €0 = lOTcSc- 



This is confirmed in Fig. 14, where the evolution of (a) the energy density, 
(b) the entropy density, (c) the pressure, and (d) the temperature is shown as a 
function of proper time t for the three equations of state {AT = 0, AT = 0.1 Tc, 
and the ideal hadron gas). Note that in the quark-gluon as well as the hadron 
phase, where p e with Cg = 1/3, Eq. (113) yields 



(118) 
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For the AT = equation of state, p = Pc = const, in the mixed phase, and (113) 
yields the cooling law 

e~r-i . (119) 

This is interpreted as follows. The longitudinal scaling expansion dilutes the sys- 
tem ^ T^^. If no mechanical work is performed, like in the mixed phase where 
dp= 0, only this geometrical dilution determines the (proper) time evolution of 
the energy density. In the phase where dp = Cg de, however, additional mechan- 
ical work is performed, and the system cools faster, e ~ r~^"'^+'^=) = t~^^^. 
The faster cooling is confirmed studying the temperature evolution, Fig. 14 
(d). For = Cg e, Cg = const., and vanishing net charges, one deduces from 
dp = c^de = sdT = {e + p)dT/T = {1 + c^)edT/T, that e ~ T^+^s"' and 
consequently, in the hadron and quark-gluon phase 

T~r-^/^ , (120) 

while in the mixed phase one deduces from dp = s dT = that 

T = const. . (121) 

This expectation is confirmed in Fig. 14 (d). 

Of course, in reality the expansion of the system will not only be purely 
longitudinal. The "Bjorkcn cylinder" will also expand transvcrsally. The principle 
of boost invariance allows us to focus on the transverse expansion at z = rj = 
only, and reconstruct the fluid properties at a different rj by performing a 
longitudinal boost with boost rapidity t/. For the sake of simplicity, let us assume 
that the system is cylindrically symmetric in the transverse direction and that 
the initial energy density profile is of the form 

e(r,ro,77 = O) = 6o0(i?,-|r|) , (122) 

where R is the transverse radius of the Bjorken cylinder. In cylindrical coor- 
dinates and a,t z = T] = 0, the conservation equations read {T^^ = E, T^^ = 
M, Vr = v): 

dtE + dr [{E + p)v] = -(^^ + l^{E+p) , (123) 

dtM + dr{Mv+p) = - (^^ + ^^ M . (124) 

Although these equations have no longer a simple analytical solution, the as- 
sumption of cylindrical symmetry has reduced the originally three-dimensional 
problem to an effectively one-dimensional problem. Indeed, for vanishing right- 
hand sides the solution of (123), (124) with the initial condition (122) is identical 
to the one of the Landau model with the substitutions z r and L ^ R. The 
right-hand sides just lead to an additional reduction of E and M from the cylin- 
drical geometry, v/r, and from longitudinal scaling, 1/t. 
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This observation, combined with the method of operator sphtting discussed 
previously, suggests the following simple solution scheme (also known as Sod's 
method [16, 29]): equations (123), (124) are of the type 

dtU + d:,F{U) = -G{U) . (125) 

The operator splitting method allows to construct the solution by first solving 
the one-dimensional partial differential equation 

dtU + d.,F{U) = Q , (126) 

(for instance with the relativistic HLLE scheme discussed above), which yields a 
prediction U for the true solution U. In a second step one corrects this prediction 
by solving the ordinary differential equation 

which is numerically realized as [30] 

U ^il - AtG{U) . (128) 

The transverse expansion of the Bjorken cylinder at z = is shown in Fig. 15 
for To — 0.1 i? and eg = 18.75 TcSc. One immediately recognizes the qualitative 
similarities with the Landau expansion, like the delay in the expansion for the 
two equations of state with a (phase) transition as compared to the expansion 
with an ideal hadron gas equation of state. The additional geometrical dilution, 
however, leads in general to a faster cooling overall and quantitatively different 
shapes for the temperature profiles and the isotherms in the t — r plane. 

Let us further quantify the time delay in the expansion induced by the tran- 
sition in the equation of state. In general, the system will decouple into free- 
streaming particles once the temperature drops below a certain "freeze-out" 
temperature Tfo, see Section 5 below. From comparison with experimental data, 
this freeze-out temperature is estimated to be on the order of 100 MeV. Let 
us therefore define a "lifetime" of the system as the time when the T = 0.7 Tc 
isotherm crosses the origin at r = in Figs. 15 (d-f). This lifetime is shown in 
Figs. 16 (a,b) as function of the initial energy density cq of the cylinder. One ob- 
serves a maximum of the such defined lifetime at initial energy densities around 
40TcSc ~ 30 GeVfm"^. At these initial energy densities, the prolongation of the 
lifetime over the respective ideal hadron gas value is about a factor of 2 (for 
AT = 0.1 Tc) to 3 (for AT = 0). 

The prolongation of the lifetime is due to the softening of the equation of 
state in the phase transition region. It is, however, interesting that the maximum 
in the lifetime does not occur around initial energy densities corresponding to 
eg (as is the case in the Landau model [23]), but at much larger initial energy 
densities. The reason for this is the strong longitudinal dilution of the system 
on account of the scaling profile Vz = z/t. In order to see a large effect of the 
softening of the equation of state in the phase transition region on the expansion 
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Fig. 15. Transverse expansion of the Bjorken cylinder for AT = (a,d), AT = 0.1 Tc 
(b,e), and the ideal gas equation of state (c,f). (a-c) show temperature profiles for 
different times, (d-f ) show the corresponding isotherms in the t — r plane (numbers are 
temperatures in units of Tc). The initial energy density is eo = 18.75 TcSc in all cases. 



dynamics, the transverse (Landau-like) expansion has to be the dominant cooling 
mechanism for the system. The Bjorken scaling expansion does not account for 
the reduced expansion tendency of the system in the transition region, it enforces 
an expansion velocity = z/t irrespective of the equation of state. In order 
to have the transverse expansion dominate the cooling of the system, one has 
to start the expansion at higher initial energy densities such that the system 
spends enough time in the mixed phase for the (slow) rarefaction shock to reach 
the origin. The initial energy density in Fig. 15 was intentionally selected to 
maximize this effect. 

Initial energy densities on the order of 10 — 30 GeVfm"^ are expected to be 
reached at the RHIC collider. In order to experimentally observe the prolongation 
of the lifetime as seen in Figs. 16, one has to find a corresponding experimental 
observable. An obvious candidate is the ratio of the "out" to the "side" radius of 
two-particle correlation functions. The "out" radius is proportional to the dura- 
tion of particle emission from a source, while the "side" radius is proportional to 
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the transverse dimension of the source (cf. [31] for a very detailed, pedagogical 
discussion). Since the transverse radius of the source is approximately the same 
in all cases, cf. Fig. 15 (a-c), the ratio i?out/-Rside seems to be a good generic 
measure for the lifctiinc. Moreover, in forming the ratio the dependence on the 
overall (unknown) spatial size of the source as well as effects from the collective 
expansion are expected to cancel. The ratio -RoutZ-Rside is plotted in Figs. 16 
(c,d) for pions with mean transverse momenta K±^ = 300 MeV. (Details on how 
to compute this quantity can be found in [30, 32].) As one observes, i?out/-Rside 
nicely reflects the excitation function of the lifetime of the system. 
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Fig. 16. Lifetime of the system as a function of eo for the Bjorken cylinder expansion, 
To = 0.1 Tc. (a) ziT = (solid) vs. ideal hadron gas (dashed), (b) AT = 0.1 Tc (dotted) 
vs. ideal hadron gas (dashed). (c,d) the corresponding ratio Rout/Rside- 
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5 Freeze-Out 

In this section I discuss an up to date unsolved problem in the application of 
relativistic fluid dynamics to describe nuclear collisions, namely the so-called 
"freeze-out" process. Given an initial condition, fluid dynamics describes the 
evolution of the system in the whole forward lightcone. Fig. 17 (a). However, 
as we have seen above, at all times near the boundary to the vacuum, as well 
as everywhere in the late stage of the evolution, the energy density becomes 
arbitrarily small, i.e., the system is rather cold and dilute. In this space-time 
region the assumption of local thermodynamical equilibrium is no longer justi- 
fied, because the particle scattering cross section a is finite, such that for small 
particle densities n the particle scattering rate, F ~ cm, becomes on the order 
of the inverse system size, F ^ . At this point, the scattering rate is too 
small to maintain local thermodynamical equilibrium and the particles decouple 
from the fiuid evolution. In this space-time region, a kinetic description for the 
particle motion would be more appropriate. One should therefore not solve fluid- 
dynamical equations in the whole forward lightcone, but only inside a space-time 
region of sufficiently large energy and particle densities, while outside this region, 
the particle motion should be described by kinetic theory. Fig. 17 (b). 



t t 




Fig. 17. (a) Conventional fluid-dynamical description in the whole forward lightcone. 
(b) Fluid dynamics describes the evolution of the system inside V4, while kinetic theory 
describes the motion of the frozen-out particles outside V4. 

The boundary S between the two regions is determined by a criterion which 
compares local scattering rates with the system size, as discussed above. The 
obvious difficulty with this more realistic description of the system's evolution is 
that this boundary has to be determined dynamically^ i.e., not only has one to 
allow for particles decoupling from the fluid, but also for the reverse process of 
particles entering the fluid from the kinetic region. (This can happen since the 
particles still, albeit rarely, collide in the kinetic region.) A consistent treatment 
of this problem is rather complicated, since one has to solve kinetic in addition 
to the fluid-dynamical equations. No serious attempt has been made so far. 
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Instead, the following approximate solution has been extensively employed: 

1. One assumes that fluid dynamics gives a reasonable description for the evo- 
lution of the system in the whole forward lightcone. 

2. One determines the "decoupling" surface E a posteriori, once the evolution 
of the fluid is known. 

3. The "thickness" of S is assumed to be infinitesimal. 

4. One assumes that particles crossing E have completely decoupled from the 
system, they stream frccily towards the detectors without any further colli- 
sional interaction ("freeze-out"). This means that they do not change their 
momentum and energy once they have crossed E. 

A very popular argument in order to determine E is the following. Since n ~ T^, 
the scattering rate F ^ (for constant cross section cr), i.e., if the temperature 
falls below a certain so-called "freeze-out" temperature Tfo, the criterion F < R 
is fulfilled, and particles decouple from the system. In this case, E is just given 
by the isotherm F = Tfo (use of this argument was already made above in the 
discussion of the "lifetime" of the system) . 

Note that assumption 3. is a strong idealization and actually rather question- 
able, because in reality E is a space-time region of finite thickness, inside which 
non-equilibrium, dissipative effects become gradually more and more important 
(the more dilute the fluid becomes), until ultimately all interactions between 
particles cease and, when leaving E, they stream freely towards the detectors. 

Nevertheless, with the above assumptions, one can readily compute the sin- 
gle inclusive spectra of particles reaching the detector. Immediately before the 
particles decouple from the fluid evolution, i.e., before they cross E, they are still 
in local thermodynamical equilibrium such that their phase space distribution 
is given by fo{k,x), Eq. (20). It is reasonable to assume that this phase space 
distribution is not changed much when they move a small distance along their 
worldlines, which carries them across E into the region of free-streaming. In that 
region, however, there are no collisions which could further change /q. There- 
fore, the phase space distribution of "frozen-out" particles is (approximately) 
the same as in local equilibrium. The total number of particles crossing a small 
surface element dE of E is then given by 

/d^k 
— di:-A:/o(fc,a;) , (129) 

N'^ being the (kinetic) particle number 4-current. The invariant momentum 
spectrum of particles crossing that surface element is consequently 

E^=dE-kfo{k,x) . (130) 

Finally, the invariant momentum spectrum (the single inclusive spectrum) of 
particles crossing the complete "freeze-out" surface E is 
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This equation is known as the Cooper-Frye formula [26], and is used in almost 
all fluid- dynamical applications to heavy-ion collisions to compute the single 
inclusive spectra of particles. 

There is, however, a problem with this formula [33]. For time-like surfaces, 
i.e., where the normal vector dSfj, is space-like, dS ■ k may either be positive 
or negative, depending on the value and direction of k^. In other words, the 
number of particles "freezing out" from a certain time-like surface clement dE 
can become negative. This is clearly unphysical, since the number of particles 
decoupling from the system must be positive definite. For space-like surfaces 
(with a time-like normal vector) as well as for time-like surface elements where 
dS ■ k > 0, the Cooper-Frye formula gives a physically reasonable, positive def- 
inite result for the number of frozen-out particles. This is illustrated in Fig. 18 
which shows the rapidity distribution of particles (i.e., the invariant momen- 
tum spectrum integrated over all transverse momenta) for massless particles 
decoupling from a freeze-out isotherm Tfo = 0.4 Tq in the Landau model with a 
p = e/3 equation of state. One clearly notices the negative particle numbers at 
midrapidity coming from the time-like parts of the isotherm. 

This contradiction is readily resolved noting that the Cooper-Frye formula 
does not really determine the number of particles decoupling from the system, 
but merely the number of particle worldlines crossing a surface element dS (and 
then integrated over the whole surface S). For time-like surface elements, there 
is of course the possibility that for certain the respective worldlines cross dS 
in the "wrong" direction, i.e., the momenta of the particles point back into the 
region of fluid, cf. Fig. 19. In particular, for the Tfo = 0.4 Tq isotherm, which 
moves away from the t-axis in the t — z plane, those are particles with vanishing 
momentum component in z direction, because their worldlines are parallel to the 
t-axis. Particles with p^ = 0, however, also have vanishing longitudinal rapidity 
y = 0, and that is the reason why these negative particle numbers appear at 
midrapidity in Fig. 18. While this explains the negative contributions in the 
Cooper-Frye formula, it also invalidates this formula as the correct prescription 
to calculate the spectra of frozen-out particles, if parts of the decoupling surface 
are time-like. 

One suggestion to circumvent this problem was to compute the final spectra 
only from contribution of particles which cross the space-like parts of E. Of 
course, as can be seen by comparing the dash-dotted with the solid line in Fig. 
18, the final spectra are dramatically different. Moreover, by neglecting particles 
crossing the time-like parts, the absolute number of frozen-out particles will also 
differ in the two cases. Note that the dN/dy distribution for particles from the 
space-like parts of the decoupling isotherm has a Gaussian shape in the Landau 
model. This was already pointed out in Landau's original paper [25] and has since 
survived as the generic (but wrong) statement that Landau's model gives rise 
to Gaussian rapidity distributions. In fact, there is no decoupling temperature 
where the full rapidity distribution including particles from the time-like parts 
resembles a Gaussian, cf. Fig. 20. 

Another suggestion to circumvent the problem of negative particle numbers 
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Fig. 18. The rapidity distribution for freeze-out along the Tfo = 0.4 To isotherm in 
the Landau model. Solid: full distribution, dotted: particles from time-like parts of the 
isotherm, dash-dotted: particles from space-like parts of the isotherm. 



is, instead of freezing out along an isotherm which has time-like parts, to freeze 
out along a surface which is space-like everywhere, for instance, a curve of con- 
stant time in the center-of-mass frame, cf. Fig. 21. In this case, all particles are 
accounted for. since the decoupling surface is bounded by the lightcone, and no 
particle can escape through the lightcone. The problem is, that also in this case, 
the spectra differ considerably from a freeze-out at constant temperature, cf. 
Fig. 22. This uncertainty is clearly unwanted when one wants to quantitatively 
compare fluid-dynamical model predictions with experimental data. 

The correct formula to compute the number of particles which physically 
decouple from the system was given in [33] : 

j dS ■kfo{k,x)0{dS -k) . (132) 

The additional 0-function ensures that negative contributions to the Cooper- 
Frye formula are cut off. The problem with this formula is that these negative 
contributions were necessary to globally conserve energy, momentum and net 
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Fig. 19. Explanation for the negative number of frozen-out particles in the Cooper- Prye 
formula. 

charge number, cf. the derivation of the conservation equations in Section 2. 
The violation of the conservation equations introduced by the freeze-out pre- 
scription (132) can, however, be circumvented by adjusting temperature, chem- 
ical potential, and the average particle 4-velocity in the single-particle distri- 
bution function fo{k,x) in (132) in such a way as to preserve the conservation 
laws. In other words, one must not use temperature, chemical potential, and 
fluid 4-velocity on the fluid side of the freeze-out surface in (132), but modi- 
fied values which ensure that energy, momentum, and net charge is conserved. 
One way to achieve this is to assume that the freeze-out surface actually is a 
conventional fluid-dynamical discontinuity across which energy, momentum, and 
net charge number are conserved. Solving the corresponding algebraic conser- 
vation equations (with energy-momentum tensor and net charge current on the 
post freeze-out side of the discontinuity constructed from (21,22) with fo{k,x) 
replaced by fo{k,x) 0{dS ■ k)) yields the required modified values for temper- 
ature, chemical potential, and average particle 4-velocity on the post freeze-out 
side. For more details, see [33, 34]. However, it still remains to be shown with 
an explicit calculation whether this suggestion to solve the freeze-out problem is 
viable in the general case. 
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